Optogenetic actuator – ERK biosensor circuits identify MAPK network nodes that shape ERK dynamics

Abstract Combining single‐cell measurements of ERK activity dynamics with perturbations provides insights into the MAPK network topology. We built circuits consisting of an optogenetic actuator to activate MAPK signaling and an ERK biosensor to measure single‐cell ERK dynamics. This allowed us to conduct RNAi screens to investigate the role of 50 MAPK proteins in ERK dynamics. We found that the MAPK network is robust against most node perturbations. We observed that the ERK‐RAF and the ERK‐RSK2‐SOS negative feedback operate simultaneously to regulate ERK dynamics. Bypassing the RSK2‐mediated feedback, either by direct optogenetic activation of RAS, or by RSK2 perturbation, sensitized ERK dynamics to further perturbations. Similarly, targeting this feedback in a human ErbB2‐dependent oncogenic signaling model increased the efficiency of a MEK inhibitor. The RSK2‐mediated feedback is thus important for the ability of the MAPK network to produce consistent ERK outputs, and its perturbation can enhance the efficiency of MAPK inhibitors.


Introduction
The extracellular signal-regulated kinase (ERK) is part of the mitogen-activated protein kinase (MAPK) signaling network and regulates a large variety of fate decisions. While ERK can be activated by several extracellular inputs, ERK signaling has mostly been studied in the context of receptor tyrosine kinases (RTKs). Upon binding of their cognate growth factors (GFs), RTKs activate a complex signaling cascade with the following hierarchy: (i) recruitment of adaptor molecules such as GRB2 (Schlessinger, 2000); (ii) control of the activity of RAS GTPases through Guanine nucleotide exchange factors (GEFs) and GTPase activating proteins (GAPs) (Cherfils & Zeghouf, 2013); (iii) triggering of a tripartite RAF, MEK, ERK kinase cascade that is further regulated by a variety of binding proteins (Lavoie et al, 2020); and (iv) ERK-mediated phosphorylation of a large number of substrates. Due to its central role in fate decisions, MAPK network dysregulation is causative for a large number of diseases including cancer (Rauen, 2013;Samatar & Poulikakos, 2014).
As for other signaling pathways (Purvis & Lahav, 2013), temporal patterns of ERK activity, hereafter referred to as ERK dynamics, rather than steady-state control fate decisions (Santos et al, 2007;Avraham & Yarden, 2011;Albeck et al, 2013;Ryu et al, 2015). These specific ERK dynamics have been shown to arise from feedback in the MAPK network. For example, a negative feedback (NFB) from ERK to RAF can produce adaptive or oscillatory ERK dynamics (Santos et al, 2007;Kholodenko et al, 2010;Avraham & Yarden, 2011). The ERK-RAF NFB was also shown to buffer against MAPK node perturbations (Sturm et al, 2010;Fritsche-Guenther et al, 2011). This property might allow cells to produce consistent ERK outputs despite heterogeneous node expression (Bl€ uthgen & Legewie, 2013). In this work, we specifically refer to the ability of the MAPK network to produce consistent ERK dynamics in presence of node perturbations as signaling robustness. While several NFBs have been mapped experimentally in the MAPK network (Lake et al, 2016), their contribution to this signaling robustness and shaping ERK dynamics remains largely unknown.
Single-cell biosensor imaging has provided new insights into MAPK signaling that were not accessible with biochemical, population-averaged measurements. This showed that the MAPK network can produce a wide variety of ERK dynamics such as transient (Ryu et al, 2015), pulsatile (Albeck et al, 2013), oscillatory (Shankaran et al, 2009), and sustained dynamics (Ryu et al, 2015;Blum et al, 2019). Mathematical modeling has provided insights into the network's structures that decode different signaling inputs into specific ERK dynamics (Santos et al, 2007;Shankaran et al, 2009;Nakakuki et al, 2010;Ryu et al, 2015). Combined modeling/ experimental approaches helped to shed light on various subparts of the MAPK network, including the epidermal growth factor receptor (EGFR) module (Koseska & Bastiaens, 2020), the RAS module (Schmick et al, 2015;Erickson et al, 2019), and the tripartite RAF/ MEK/ERK cascade (Ferrell & Bhatt, 1997;Kholodenko, 2000;Orton et al, 2005;Santos et al, 2007;Ryu et al, 2015;Kochańczyk et al, 2017;Arkun & Yasemi, 2018). However, the low experimental throughput to measure ERK dynamics, or other MAPK network nodes, has precluded a global understanding of the specific functions of the nodes present in the network.
Here, we built multiple genetic circuits consisting of optogenetic actuators together with an ERK biosensor to simultaneously activate ERK from different nodes in the MAPK network and report singlecell ERK dynamics. These circuits allowed us to investigate the role of 50 MAPK signaling nodes in ERK dynamics regulations with RNA interference (RNAi). We observed that most perturbations of individual nodes resulted in mild ERK dynamics phenotypes despite targeting major MAPK signaling nodes. Further, the ERK dynamics induced by various perturbations suggest that two NFBs (ERK-RAF and ERK-RSK2-SOS) act simultaneously to regulate ERK dynamics. Targeting the RSK2-mediated NFB increased the efficiency of additional MAPK network perturbations both in our optogenetic systems and in an ErbB2-driven oncogenic ERK signaling model. This suggests that the RSK2-mediated feedback plays a role in MAPK signaling robustness and can be targeted for potent inhibition of oncogenic ERK signaling.

Results
An optogenetic actuator-biosensor genetic circuit to study input-dependent ERK dynamics In order to measure ERK dynamics in response to dynamic RTK input, we built a genetically encoded circuit made of an optogenetic RTK actuator and an ERK biosensor ( Fig 1A). We chose optoFGFR, which consists of a myristoylated intracellular domain of the fibroblast growth factor receptor 1 (FGFR1) fused to a CRY2 domain and tagged with mCitrine . Upon stimulation with blue light, optoFGFR dimerizes and trans-autophosphorylates, leading to the activation of the MAPK/ERK, phosphoinositide 3-kinase (PI3K)/ AKT, and phospholipase C (PLC)/Ca 2+ pathways. As ERK biosensor, we used ERK-KTR-mRuby2 that is spectrally compatible with optoFGFR. ERK-KTR reversibly translocates from the nucleus to the cytosol upon ERK activation (Regot et al, 2014). We used a nuclear Histone 2B (H2B)-miRFP703 marker to identify and track single cells. After stably inserting these constructs into murine NIH3T3 fibroblasts, we used automated time-lapse microscopy to stimulate selected fields of view with defined blue light input patterns to activate optoFGFR. The corresponding ERK-KTR/H2B signals were recorded with a 1-min temporal resolution. We observed that a 100ms light pulse leads to reversible ERK-KTR translocation from the nucleus to the cytosol, indicative of transient ERK activation (Fig 1  B, Movie EV1). At the end of each experiment, we imaged the mCitrine signal to evaluate optoFGFR expression levels. We built a computer vision pipeline to automatically track each nucleus, compute ERK activity as the cytosolic/nuclear ratio of the ERK-KTR signals, and correlate single-cell ERK responses with optoFGFR levels (Fig 1C). We then use this pipeline to evaluate the sensitivity and specificity of our system with dose-response experiments using the FGFR inhibitor SU5402, the RAF inhibitor RAF709, the MEK inhibitor U0126, and the ERK inhibitor SCH772984 (Appendix Fig S1A).
To evaluate light-dependent optoFGFR activation dynamics, we engineered a mScarlet-tagged optoFGFR that is spectrally orthogonal to CRY2 absorption (Appendix Fig S1B). Total internal reflection (TIRF) microscopy visualized the formation of optoFGFR clusters in response to blue light-mediated dimerization in the plasma membrane (Appendix Fig S1B, blue arrows, Movie EV2). Consistently with CRY2's dissociation half-life (Duan et al, 2017), these optoFGFR clusters appeared within 20 s after a blue light pulse and disappeared after~5 min (Appendix Fig S1C). We assume that optoFGFR is active in its clustered form in which transphosphorylation occurs and inactive in its monomeric form due to tonic cytosolic phosphatase activity (Lemmon et al, 2016). As documented previously , light stimulation also triggered optoFGFR endocytosis (Appendix Fig S1B, red arrows).
Directly following light stimulation, we systematically observed a short ERK inactivation period, that we refer to as "dip", lasting 2-3 min before activation of a strong ERK activity (Appendix Fig S1D, green rectangle). This light-induced ERK dip was insensitive to SCH772984-mediated ERK inhibition but could be suppressed by cyclosporine A-mediated calcineurin inhibition. Calcineurin is a Ca 2+dependent phosphatase that dephosphorylates Ser383 in Elk1 (Sugimoto et al, 1997). As ERK-KTR contains an Elk-1 docking domain phosphorylated by ERK (Regot et al, 2014), we hypothesized that it could be negatively affected by optoFGFR-evoked Ca 2+ input ) (Appendix Fig S1E). Consistently, Ionomycin-evoked increase in cytosolic Ca 2+ induced a dip in absence of light stimulation (Appendix Fig S1F).

Different optoFGFR inputs trigger transient, oscillatory, and sustained ERK dynamics
Next, we characterized optoFGFR-triggered ERK dynamics in response to a single light pulse of different intensities and durations (Fig 2A). As ERK dynamics depended on light power density, as well as pulse duration, we defined the light dose (D, mJ/cm 2 ) as their product to quantify the total energy received per illuminated area. To characterize ERK dynamics, we extracted the amplitude at the maximum of the peak (maxPeak), and the full width at half maximum (FWHM) of the ERK trajectories ( Fig 2B). With increasing light doses, ERK peaks increased both in duration and amplitude, until the latter reached saturation. Based on these observations, we selected 180 mW/cm 2 and 100 ms (D = 18 mJ/cm 2 ) as the minimal light input to generate an ERK transient of maximal amplitude. Using this light dose, we then investigated ERK dynamics in response to multiple light pulses delivered at different intervals ( Fig 2C). All stimulation regimes led to identical maximal ERK amplitude ( Fig EV1A) and adaptation kinetics when optoFGFR input ceased ( Fig EV1B). Repeated light inputs applied at 10-or 20-min intervals evoked population-synchronous ERK transients. By contrast, repeated light inputs applied at higher frequencies (2-min intervals) led to sustained ERK dynamics. Given CRY2's 5-min dissociation half-life (Appendix Fig S1B and C) (Duan et al, 2017), this suggests that light pulses delivered at a 2-min interval reactivate optoFGFR faster than it deactivates, leading to sustained optoFGFR activity. Hierarchical clustering of ERK responses to sustained optoFGFR input highlighted the presence of sustained and oscillatory single-cell ERK dynamics (Fig 2D). Classification of ERK trajectories based on optoFGFR expression revealed that sustained/ oscillatory ERK dynamics correlated with high/low optoFGFR levels ( Fig 2E, Movie EV3). Oscillatory ERK dynamics were also observed in optoFGFR high expressing cells in response to low light input ( Fig 2F). Thus, sustained optoFGFR input can trigger sustained or oscillatory ERK dynamics depending on the input strength, a combination of light energy and optoFGFR expression.

ERK dynamics evoked by optoFGFR versus endogenous RTKs highlight different MAPK regulatory mechanisms
Because of the absence of an ectodomain, optoFGFR must be considered as a prototypic RTK that lacks some regulatory mechanisms inherent to the native FGFR. To evaluate whether optoFGFR is relevant for studying the MAPK network, we compared ERK dynamics evoked by optoFGFR inputs versus stimulation of the endogenous FGFR or EGFR using increasing concentrations of basic FGF (bFGF) and EGF. All bFGF concentrations led to an ERK peak similar in amplitude to sustained optoFGFR input (Figs 3A and EV1C, compared to EV1A). However, FGFR inputs led to different ERK dynamics than optoFGFR: 1 ng/ml bFGF led to damped ERK oscillations followed by steady-state sustained ERK activity, while 10 and 100 ng/ml bFGF concentrations led to a first ERK peak followed by a strong adaptation. The biphasic behavior induced by increasing bFGF concentrations was previously documented to emerge from the competition of bFGF for FGFR and heparan sulfate proteoglycan co-receptors (Kanodia et al, 2014;Blum et al, 2019). It is thus not surprising that optoFGFR, that lacks these extracellular interactions, produced different ERK dynamics than FGFR. All EGF concentrations led to an ERK peak similar in amplitude to optoFGFR and FGFR inputs (Figs 3B and EV1D Figure 1. An optogenetic actuator-biosensor genetic circuit to study input-dependent ERK dynamics. A Schematic representation of the optoFGFR system consisting of the optogenetic FGF receptor (optoFGFR) tagged with mCitrine, the ERK biosensor (ERK-KTR) tagged with mRuby2, and a nuclear marker (H2B) tagged with miRFP703. B Time-lapse micrographs of ERK-KTR dynamics in response to a 470-nm light pulse. Using a 20× air objective, ERK-KTR and H2B channels were acquired every 1 min and the optoFGFR channel was acquired once at the end of the experiment. Scale bar: 50 μm. C Image analysis pipeline developed to quantify single-cell ERK dynamics. Nuclear and cytosolic ERK-KTR signals were segmented based on the H2B nuclear mask.
Single-cell ERK activity was then calculated as the cytosolic/nuclear ERK-KTR ratio. Single-cell optoFGFR intensity was measured under the cytosolic ERK-KTR mask and used as a proxy for single-cell optoFGFR expression.   To provide intuition about the MAPK network circuitries leading to different ERK dynamics in response to optoFGFR and EGFR inputs, as well as the origin of the oscillatory behavior, we built a mathematical model consisting of the RAS GTPase and the threetiered RAF/MEK/ERK network ( Fig 3D, Appendix Table S1). We used ordinary differential equations with Michaelis-Menten kinetics (see Materials and Methods, Appendix Table S2 and S3). To account for the oscillatory ERK dynamics in response to EGFR and optoFGFR inputs, we included the well-documented ERK-RAF NFB (Santos et al, 2007;Kholodenko et al, 2010;Fritsche-Guenther et al, 2011;Blum et al, 2019). We also included a receptor-level inactivation process for EGFR, but not for optoFGFR, to account for EGFdependent regulatory mechanisms. We used a Bayesian inference approach (Mikelson & Khammash, 2020) to infer the model parameters from averaged ERK trajectories in response to sustained low optoFGFR input with or without sustained EGFR pre-stimulation ( Fig 3E). After identification of parameters that allowed the model to capture the training dataset ( Fig 3F), we simulated ERK dynamics evoked by low EGFR input (adaptative, oscillatory ERK dynamics), high EGFR input (adaptative ERK dynamics without oscillation), and sustained high optoFGFR input (sustained ERK dynamics) ( Fig 3G). We observed that our model with a NFB and EGFR inactivation was able to predict ERK dynamics evoked by different EGFR and optoFGFR input strengths, while two simpler models (one with only the EGFR inactivation reaction, but no NFB (Fig EV1E-G) and one with only the NFB, but no EGFR inactivation (Fig EV1H-J) were not able to reproduce experimentally observed ERK dynamics.
This suggested that oscillatory optoFGFR-evoked ERK dynamics emerge from a NFB also present downstream of endogenous EGFR, while additional regulatory mechanisms seem to be required for the strong ERK transient adaptation following EGFR input. These mechanisms might consist of receptor-level regulations such as endocytosis, which was recently shown to be an important regulator of the transient adaptive EGF-triggered ERK dynamics in different cell systems (Gerosa et al, 2020;Kiyatkin et al, 2020). While optoFGFR also gets endocytosed (Appendix Fig S1B, Kim et al, 2014), it most likely is insensitive to inactivation by endosome acidification since it lacks an ectodomain (Huotari & Helenius, 2011). Additionally, light-mediated optoFGFR dimerization might occur both at the plasma and endo-membranes, allowing for reactivation of endocytosed optoFGFR. The hypothesis that a receptor-level mechanism is important for strong adaptation was further supported by inhibition of optoFGFR with the FGFR kinase inhibitor (SU5402), which shifted ERK dynamics from sustained to transient in a dose response-dependent manner ( Fig EV1K). Thus, these results suggest that optoFGFR lacks receptor-dependent regulatory mechanisms but allows us to investigate the intracellular MAPK feedback structure shaping ERK dynamics. In our model, we used the well-established ERK-RAF NFB. However, several NFBs have been mapped in the MAPK signaling cascade, whose role in shaping ERK dynamics is still unknown and which could also be responsible for the observed oscillatory ERK dynamics.
RNA interference screen reveals that ERK dynamics remain unaffected in response to perturbation of most MAPK signaling nodes We then explored the network circuitry that shapes optoFGFRevoked ERK dynamics with an RNA interference (RNAi) screen targeting 50 MAPK signaling nodes. We focused our screen on sustained optoFGFR input which captured the largest amount of information about ERK dynamics when compared to other stimulation schemes: it led to sustained and oscillatory ERK dynamics (Fig 2E and F) while recapitulating the rapid increase of ERK activity and adaptation observed with transient input (Fig EV1A and B). We used a bioinformatic approach to select 50 known interactors of the tripartite RAF/MEK/ERK cascade downstream of the FGFR receptor that were detected in a NIH3T3 proteome (Schwanh€ ausser et al, Table S4). We used the siPOOL technology to specifically knockdown (KD) these 50 MAPK signaling nodes while limiting off-target effects (Hannus et al, 2014). We first validated KD efficiency by quantifying transcript levels with different siPOOL concentrations targeting the ERK and MEK isoforms ( Fig  EV2A) and observed strong KD with 10 nM siRNA concentration. We then evaluated the effect of ERK1 or ERK2 KD on ERK dynamics. We observed only subtle phenotypes compared to the non-targeting siRNA (CTRL) used as negative control (Fig 4B), even though efficient KD was observed at protein level ( Fig 4C). However, combined ERK1/ERK2 KD strongly suppressed ERK dynamics indicating that the latter is not affected by the perturbation of individual ERK isoforms as previously reported (Fritsche-Guenther et al, 2011;Ornitz & Itoh, 2015). Due to its strong phenotype, we used ERK1/ ERK2 KD as positive control throughout our screen. A ERK responses to increasing light power densities and pulse durations of 470-nm transient light input. The light dose "D" is calculated as the product of the power density and pulse duration. B Quantification of the maxPeak (maximal ERK amplitude of the trajectory) and the FWHM (full width at half maximum) of single-cell ERK responses shown in (A) (N min = 40 cells per condition from one replicate, boxes indicate the upper and lower quartiles, the central bands indicate the median and whiskers extend to individuals up to 1.5 interquartile away from the median). C ERK responses to 470 nm light pulses delivered every 20, 10, 5, and 2 min, respectively (D = 18 mJ/cm 2 ). D Hierarchical clustering (Euclidean distance and Ward D2 linkage) of trajectories from the 2-min interval stimulation shown in (C) (referred to as "sustained") (N = 60 cells). The number of clusters was empirically defined to resolve the different ERK dynamics. The average ERK responses per cluster are displayed on the right. E Using the log10 intensity of optoFGFR-mCitrine (left panel), trajectories shown in (D) are separated in low and high optoFGFR cells (right panel). F ERK responses to increasing doses of sustained optoFGFR input. Single-cell ERK trajectories were divided in low (top panel) and high (bottom panel) optoFGFR expression.
ª 2022 The Authors Molecular Systems Biology 18: e10670 | 2022 We performed three replicates of the screen targeting the 50 nodes. Despite efficient KD quantified for different nodes (Fig EV2B), visual inspection of ERK trajectories only revealed subtle ERK dynamics phenotypes for a limited number of node perturbations ( Fig EV2C and D). We used a feature-based approach to evaluate the effect of each perturbation on ERK dynamics. We focused our analysis on ERK responses evoked by high optoFGFR input to limit the single-cell heterogeneity due to optoFGFR expression variability.    High optoFGFR We quantified the average ERK activity before stimulation (baseline), the maximal ERK amplitude during stimulation (maxPeak), and the ERK amplitude at a fixed time point after response adaptation in the negative control (ERKpostStim). To evaluate these phenotypes, we z-scored the features associated to each perturbation to those of the negative control ( Fig 4D, see Materials and Methods for details). While many phenotypes were statistically significant, most of them remained mild as observed by visually inspection of the feature distributions ( Fig EV3A). Apart from ERK1+2 KD, only GRB2, PTK2, and ERK2 led to a reduction of ERK amplitude (maxPeak). KD of negative regulators such as SPROUTY 2,3, and 4, or phosphatases such as PP2A and several dual-specificity phosphatases (DUSPs) led to increased ERK amplitude. Increased basal ERK activity was observed for RKIP, PP2A, DUSP4, and DUSP6 KDs, indicating a function in regulating basal ERK levels. Prolonged ERK activity (ERKpostStim) was observed in response to KD of RKIP, PP2A, ERK2, DUSP1,2,3,4,6 and strikingly for RSK2 KD (Fig EV3B), suggesting a role of these nodes in ERK adaptation. Because both visual inspection of trajectories, as well as our feature-based approach might miss more subtle ERK dynamics phenotypes, we used CODEX (Jacques et al, 2021), a data-driven approach to identify patterns in single-cell time series based on convolutional neural networks (CNNs) (Fig EV3C). We trained a CNN to classify ERK trajectories that originate from different siRNA perturbations and selected the ten perturbations for which the CNN classification accuracy was the highest (Appendix Table  S4, "CODEX accuracy", see Materials and Methods for details). Projection of the CNN features in a t-distributed stochastic neighbor embedding (t-SNE) space revealed different clusters of ERK trajectories ( Fig EV3D). Comparison of the ten trajectories with the highest classification confidence identified by CODEX to randomly selected ERK trajectories for low or high optoFGFR expression highlighted ERK phenotypes not accessible to visual inspection and the feature-based approach (Fig 4E). CODEX identified some of the perturbations that affect ERK amplitude, baseline or adaptation observed with the feature-based approach. However, it also highlighted perturbations affecting oscillatory ERK dynamics. PP2A KD led to sustained oscillatory behavior. PLCG1 KD resulted in a first peak followed by damped oscillations, and absence of the dip. As phospholipase C mediates Ca 2+ signaling in response to FGFR activation (Ornitz & Itoh, 2015), this further validates the role of Ca 2+ signaling in formation of the dip (Appendix Fig S1D-F). RAPGEF1 KD led to oscillatory ERK responses of different amplitudes. RSK2, ERK2, and CRAF KD displayed reduced oscillatory ERK behavior.
To validate the latter oscillatory ERK dynamics phenotypes, we evaluated the proportion of oscillatory trajectories (trajectories with at least 3 peaks) for each perturbation, both for high and low optoFGFR input ( Fig 4F). This confirmed that RSK2, CRAF, and ERK2 KD led to decreased oscillatory ERK dynamics. We also observed that these perturbations reduced ERK oscillations in cells stimulated with 1 ng/ml EGF (Fig EV3E-G), suggesting a role of these nodes in the regulation of ERK oscillations in the context of a native RTK.
ERK2 and CRAF isoforms are implicated in the well-established ERK-RAF NFB, known to regulate ERK dynamics (Santos et al, 2007;Ryu et al, 2015;Blum et al, 2019), and to enable consistent ERK dynamics under MEK or ERK perturbations (Sturm et al, 2010;Fritsche-Guenther et al, 2011). RSK2 encodes the p90 ribosomal S6 kinase 2 protein, an ERK substrate regulating survival and proliferation (Cargnello & Roux, 2011;Yoo et al, 2015). RSK2 is also known to be involved in an ERK-induced NFB targeting SOS (Douville & Downward, 1997;Saha et al, 2012;Lake et al, 2016), whose significance in the regulation of ERK dynamics has been less well studied. In addition to dampening ERK oscillations, RSK2 KD also led to slower ERK adaptation when optoFGFR input ceased (Figs 4D and EV3A and B), suggesting an important role of this NFB in ERK dynamics regulation. Our results suggest that the ERK-RAF and ERK-RSK2-SOS NFBs simultaneously operate within the MAPK network to generate ERK oscillations and raise the question whether both NFBs contribute to the strong MAPK signaling robustness observed in our screen.

Direct optogenetic activation of RAS highlights different ERK dynamics phenotypes than optoFGFR input
To further explore the role of MAPK feedback in MAPK signaling robustness, we used optoSOS (Johnson et al, 2017), an optogenetic actuator that activates RAS, and thus bypasses the RSK2-mediated NFB regulation (Fig 5A). OptoSOS consists of a membrane anchored light-activatable iLID domain, and an mCitrine-tagged SspB domain fused to SOS's catalytic GEF domain. It was stably integrated into ◀ Figure 4. RNA interference screen reveals that ERK dynamics remain unaffected in response to perturbation of most MAPK signaling nodes.
A RNAi perturbation targets referred to by their protein names. Nodes were spatially grouped based on the hierarchy of interactions within the MAPK network and color-coded for their function. B ERK responses to sustained optoFGFR input (D = 18 mJ/cm 2 ) in cells transfected with 10 nM siRNA against ERK1, ERK2 or a 5 nM combination of each (ERK1+2). A non-targeting siRNA (CTRL) was used as control (N = 15 cells from low and high optoFGFR levels). C Western blot analysis of cells transfected with 10 nM siRNA against ERK1, ERK2 or a 5 nM combination of each (ERK1+2). D Z-Score evaluation of the baseline, maxPeak and ERKpostStim of single-cell ERK responses under sustained high optoFGFR input (D = 18 mJ/cm 2 ). The z-score was calculated by comparing each RNAi perturbation to the CTRL KD (N min = 126 cells per treatment, from three technical replicates). Non-significant (NS) results are in gray (see Fig EV3A for statistical results). E Single-cell ERK trajectories (sustained optoFGFR input, D = 18 mJ/cm 2 ) for the RNAi perturbations classified with the highest accuracy by CODEX. Top lines show single-cell ERK trajectories for which CODEX had the highest classification confidence in the validation set (N = 10). Bottom lines show single-cell ERK trajectories for low and high optoFGFR cells (N = 30 for each condition, randomly selected out of at least 212 cells per perturbation from three technical replicates). For easier visualization, the CTRL condition is shown twice. F Proportion of oscillating cells (trajectories with at least three peaks) per RNAi perturbation for low and high optoFGFR expression (sustained optoFGFR input, D = 18 mJ/cm 2 , N min = 61 cells for low and 126 for high optoFGFR per perturbation from three technical replicates). Perturbations were ordered based on the proportion of oscillating cells with low optoFGFR expression. Statistical analysis was done using a pairwise t-test, comparing each perturbation against the CTRL for each receptor level independently (*< 0.05, **< 0.005, ***< 0.0005, ****< 0.00005, FDR P-value correction method).
cells expressing ERK-KTR and H2B. Because iLID displays faster dissociation rates than CRY2 (t 1/2 = 30 s for iLID versus~5 min for CRY2 (Duan et al, 2017;Benedetti et al, 2018)), optoSOS required repeated light pulses to prolong its membrane recruitment and produce a robust ERK response ( Fig 5B). Five consecutive 100 ms light pulses at 6 W/cm 2 (D = 0.6 J/cm 2 ) applied at 20-s intervals, provided the minimal light input to induce a saturated ERK amplitude ( Fig EV4A). Application of this light input at 2-min intervals evoked sustained ERK dynamics with small fluctuations at the same frequency as the light input pattern, reflecting the fast optoSOS  reversion to the dark state ( Fig 5C). OptoSOS did not induce ERK oscillations ( Fig EV4B), even in cells with low optoSOS expression or at lower light doses ( Fig 5D). However, ERK amplitudes correlated with optoSOS expression level, low optoSOS levels led to low ERK amplitudes, while high actuator expression levels resulted in high ERK amplitudes. Using the minimal light input to trigger saturating ERK amplitude, both optoSOS and optoFGFR led to steep ERK activation and fast adaptation when light stimulation ceased (compare Figs 2C and 5C), as well as similar ERK amplitudes in cells expressing high actuator levels ( Fig 5E). However, high optoSOS expression levels moderately increased ERK activity baseline levels in comparison to optoFGFR (Fig EV4C), suggesting that this system is leaky to some extent. Using this specific light input, we performed siRNA screens targeting MAPK signaling nodes downstream of optoSOS in triplicates (Fig EV4D and E). We extracted the baseline, maxPeak, ERKpostStim features from optoSOS high expressing cells ( Fig  EV4F) and z-scored feature values to the negative control ( Fig 5F). We observed more prominent ERK amplitude phenotypes in response to optoSOS input than to optoFGFR input. Some of these phenotypes are shown in Fig 5G. Most prominently, CRAF, ERK2, and DUSP4 KD led to a stronger reduction in ERK amplitude than observed with optoFGFR input. RSK2 KD also reduced ERK amplitude, suggesting that it also regulates nodes downstream of RAS. However, RSK2 KD did not decrease ERK adaptation following optoSOS input removal (Fig EV4G), suggesting that it is not involved in NFB regulation in this system. PP2A KD did not induce increased ERK amplitude or baseline as observed in the optoFGFR system. As for optoFGFR input, DUSP6 KD increased basal ERK activity and decreased adaptation (Fig EV4G). DUSP22 KD led to increased amplitude, without affecting ERK baseline and adaptation. NF1 KD, which encodes a RAS-specific GAP, led to increased ERK baseline and slower adaptation (Fig EV4G), without affecting ERK amplitude. The NF1 baseline phenotype, that was not observed in the optoFGFR system, might emerge from the optoSOS-mediated low levels of RAS activation due to the optoSOS system's leakiness (Fig EV4C) that can then be amplified by loss of NF1's RAS GAP activity. The finding that perturbation of specific nodes (e.g., ERK2 and CRAF) leads to more penetrant phenotypes in response to optoSOS versus optoFGFR input suggested that the RAS/RAF/MEK/ERK part of the network is more sensitive to perturbations than optoFGFR-triggered network, suggesting that the RSK2 NFB that operates above RAS contributes to MAPK signaling robustness.

Perturbation of the RSK2-mediated NFB increases the efficiency of RAF, MEK, and ERK targeting drugs
To further investigate the role of the RSK2-mediated NFB in MAPK signaling robustness, we performed dose response experiments using different MAPK inhibitors and compared ERK amplitudes evoked by optoFGFR (RSK2-feedback dependent) versus optoSOS (RSK2-feedback independent) input, as well as optoFGFR input in absence/presence of RSK2 perturbation. We used drugs targeting B/ CRAF (RAF709), MEK (U0126) and ERK (SCH772984). We evaluated the inhibition efficiency by measuring ERK amplitude at a fixed time point, focusing on ERK responses evoked by high optoFGFR or optoSOS inputs to limit the single-cell heterogeneity due to expression variability of the optogenetic actuator. All inhibitors led to a stronger reduction of ERK amplitude and EC 50 in response to optoSOS versus optoFGFR input (Figs 6A-C and EV5A, Appendix Table S5). Visual evaluation of ERK amplitude distributions ( Fig 6B) and quantification of their standard deviations ( Fig 6D) revealed more compact ERK amplitude distributions in presence of increasing drug concentrations in response to optoSOS versus optoFGFR input. This suggests a more homogeneous drug inhibition in the cell population in response to optoSOS input. We then performed the identical experiments in CTRL or RSK2 KD cells in response to optoFGFR input (Figs 6E-H and EV5B, Appendix Table S6). RSK2 KD led to increased inhibition of ERK amplitudes, decreased EC 50 , and more compact ERK amplitude distributions in response to increasing drug concentration than in CTRL KD cells. Similar results were observed when the RSK2-mediated feedback was inhibited using the RSK inhibitor SL0101 (Smith et al, 2005) (Fig EV5C-F, Appendix Table S7). Thus, inhibition of the RSK2-mediated NFB sensitizes ERK responses to RAF, MEK or ERK drug perturbations. Note that drug mediated ERK amplitude inhibition was stronger in response to optoSOS input than to optoFGFR input with RSK2 KD or RSK inhibition, suggesting that additional mechanisms to the RSK2mediated feedback contribute to MAPK signaling robustness. However, our results suggest that perturbation of the RSK2-mediated feedback can be exploited to enhance the efficiency of MAPK-targeting drugs, reducing ERK amplitudes more homogeneously across the cell population.
◀ Figure 5. Direct optogenetic activation of RAS highlights different ERK dynamics phenotypes than optoFGFR input.
A Schematic representation of ERK signaling induced by optoSOS versus optoFGFR input. B ERK dose responses under transient optoSOS input consisting of different numbers of repeated 470 nm pulses (1, 2, 3, 4, and 5× pulses applied at 20-s intervals, D = 0.6 J/cm 2 ). Repeated pulses are depicted as a single stimulation (blue bar). C ERK responses to optoSOS inputs consisting of 5 repeated 470 nm light pulses delivered every 20, 10, 5, and 2 min, respectively (D = 0.6 J/cm 2 ). D ERK responses to increasing light doses of sustained optoSOS input consisting of 2-min interval input, each input made of 5 repeated light pulses. Cells were divided in low and high optoSOS expression levels based on the log10 intensity of the optoSOS-mCitrine. E Quantification of the maxPeak of single-cell ERK responses under sustained optoFGFR (Fig 2F, D = 18 mJ/cm 2 ) and optoSOS (Fig 5D, D = 0.6 J/cm 2 ) input for low or high expression of each optogenetic system (N = 40 cells per condition, boxes indicate the upper and lower quartiles, the central bands indicate the median and whiskers extend to individuals up to 1.5 interquartile away from the median). Statistical analysis was done using a Wilcoxon test, comparing each condition to each other (N min = 48 cells per condition, NS: non-significant, *< 0.05, ****< 0.00005, FDR P-value correction method). F Z-Score evaluation of the baseline, maxPeak, and ERKpostStim of single-cell ERK responses under sustained high optoSOS input (D = 0.6 J/cm 2 ). The z-score was calculated by comparing each RNAi perturbation to the CTRL KD (N min = 33 cells per treatment, from three technical replicates). Non-significant (NS) results are in gray (see Fig EV4F for  Targeting the RSK2-mediated feedback in an ErbB2 oncogenic signaling model increases MEK inhibition efficiency The results above suggested an important role of the RSK2-mediated feedback in MAPK signaling robustness against node perturbation in response to optogenetic inputs in NIH3T3 cells. To test whether this feedback also contributes to MAPK signaling robustness in a disease-relevant system, we evaluated its function in MCF10A cells, a breast epithelium model, using either wild-type (WT) or overexpressing ErbB2 (referred to as ErbB2 over ) recapitulating the ErbB2 amplification observed in 20% of all breast cancers (Yarden & Pines, 2012;Arteaga & Engelman, 2014). We chose this specific model system because ErbB2 amplification leads to constitutive RTK input on the MAPK network, while retaining an intact downstream feedback structure (Fig 7A). This contrasts with other cancer model systems in which additional mutations might lead to RAS or RAF overactivation, and thus disrupt the feedback architecture. Further, previous work has highlighted the role of NFBs in ERK pulse formation in MCF10A cells (Kochańczyk et al, 2017), suggesting that EGFR and ErbB2 trigger a MAPK network with similar feedback circuitry as optoFGFR.
As described before (Albeck et al, 2013), WT cells displayed asynchronous low-frequency ERK pulses in the absence of EGF, and high-frequency ERK pulses in presence of EGF (Fig 7B). In marked contrast, ErbB2 over cells displayed high frequency ERK pulses, even in the absence of EGF (Fig 7C). To investigate the role of the RSK2mediated feedback in MAPK signaling robustness, we performed a U0126 dose response in EGF-stimulated MCF10A WT cells and found that 3 µM U0126 decreased ERK amplitude without fully suppressing the response (Fig EV5G and H). As observed in response to optogenetic inputs, RSK inhibition with 50 µM SL0101 led to a mild reduction in ERK amplitude. However, in combination with 3 µM U0126, ERK amplitude was decreased to the level of unstimulated cells. Similar results were observed in ErbB2 over cells (Fig EV5I), suggesting that RSK2 perturbation increases the sensitivity of ERK responses to MEK inhibition.
As averaging ERK dynamics can hide asynchronous single-cell signaling activity, we further investigated the effect of these perturbations on single-cell trajectories using CODEX (Jacques et al, 2021) (see Materials and Methods for details). For WT cells, a tSNE projection of the CNN features built from single-cell ERK trajectories hinted that the CNN was able to construct features separating the treatments into well-defined clusters (Figs 7D and EV5J). Clustering of the CNN features confirmed the existence of discrete ERK dynamics clusters (Fig 7E) whose composition correlated with the treatments (Fig 7F). To characterize the dynamics captured by each cluster, we extracted the medoid trajectory and its 4 closest neighbors from each cluster (Fig 7G). This revealed that non-stimulated cells mostly display low frequency ERK activity pulses (cluster 4) or absence of pulses (cluster 5). Cells stimulated with EGF without inhibitor displayed ERK pulses of high amplitude (cluster 1). SL0101-treated cells displayed a sustained ERK activation at low amplitude (cluster 3). U0126-treated cells still displayed prominent ERK pulses but at a lower amplitude than EGF-treated cells in absence of drug (cluster 2). Finally, in cells treated with both U0126 and SL0101, almost no ERK activity was observed (cluster 5). For ErbB2 over cells, we observed that the CNN features were forming a more continuous space with less distinct clusters (Figs 7H and I, and EV5K). A heterogeneous mix of ERK trajectory clusters was observed for the different treatments (Fig 7J and K). Untreated cells mostly displayed high frequency ERK pulses that were either sharp (cluster 3) or wider (cluster 2). SL0101-treated cells were almost equally shared between cluster 1 (relatively flat high-amplitude ERK trajectories), cluster 2, cluster 4 (low-amplitude ERK pulses) and cluster 5 (low baseline activity). U0126 led to a less heterogeneous mix mostly consisting of ERK trajectories from cluster 4 and 5. The U0126/SL0101 combination shifted most cells to cluster 5, indicating an efficient inhibition of ERK activity at a suboptimal U0126 concentration. Thus, RSK inhibition also sensitizes the MAPK network to U0126-mediated MEK inhibition both in MCF10A WT and ErbB2over cells.

Discussion
Optogenetic actuator-biosensor circuits allow for feedback structure mapping in the MAPK network ERK dynamics is crucial for fate decisions. Yet, the topology of the network enabling the cells to sense different inputs and convert this information into finely tuned ERK dynamics remains poorly understood. We developed genetic circuits consisting of optogenetic actuators and an ERK biosensor (Figs 1A and 5A) that allow for a largescale interrogation of single-cell ERK dynamics and investigated the ▸ Figure 7. Targeting the RSK2-mediated feedback in an ErbB2 oncogenic signaling model increases MEK inhibition efficiency.

A
Schematic representation of MAPK signaling in response to EGFR input in MCF10A WT cells or oncogenic ErbB2 input in ErbB2 overexpressing (ErbB2 over ) cells. B, C Single-cell ERK responses in MCF10A WT cells without or with stimulation with 10 ng/ml EGF at t = 30 min (B) and in unstimulated MCF10A ErbB2 over cells (C). D tSNE projection of CODEX's CNN features from ERK trajectories of MCF10A WT cells without EGF stimulation, or with 10 ng/ml EGF stimulation added at t = 30 min in absence of perturbation, with 50 μM SL0101, 3 μM U0126 or a combination of both. effects of 50 RNAi perturbations targeting components of the MAPK signaling network (Fig 4A). In our optoFGFR screen, we only observed a small number of penetrant ERK dynamics phenotypes (Fig 4D-F), implying that the MAPK network can buffer against perturbations of most of its components. We cannot exclude that in some cases, even on the relatively short 72 h timescale of the RNAi experiment, compensation by upregulation of specific nodes might occur. However, our data suggest that the MAPK network topology allows for MAPK signaling robustness-the production of consistent ERK outputs in presence of node perturbation. This might emerge from isoform redundancy for multiple nodes in the network, as observed for single or combined ERK isoforms perturbation (Fig 4B)  decrease the network sensitivity to node perturbation (Sturm et al, 2010;Fritsche-Guenther et al, 2011). Our screen suggested that RSK2, that mediates a NFB from ERK to SOS (Douville & Downward, 1997;Saha et al, 2012), both regulates ERK dynamics (Fig 4D-F) and plays a role in MAPK signaling robustness (Fig 6E-H).
In addition, our data suggest that the well-studied ERK-RAF NFB, which has been shown to buffer against MAPK node perturbations (Sturm et al, 2010;Fritsche-Guenther et al, 2011), also regulates ERK dynamics (Fig 4F). We therefore speculate that feedback operates simultaneously in the MAPK network, and act at multiple levels within the cascade to warrant MAPK signaling robustness. Consistently with this hypothesis, we observed that the optoSOStriggered network, which is not under the RSK2 NFB regulation, shows an increased sensitivity in ERK amplitude to perturbation of some nodes (Fig 5F and G). Indeed, ERK2 and CRAF perturbations, which led to loss of ERK oscillations, had relatively mild amplitude phenotypes in response to optoFGFR input, while both perturbations led to strong ERK amplitude phenotypes in response to optoSOS input. Because these phenotypes were not observed with other ERK and RAF isoforms, we propose that ERK2 and CRAF are the isoforms involved in the classic ERK-RAF NFB. Additional feedback has been reported within the MAPK network (Langlois et al, 1995;Lake et al, 2016;Kochańczyk et al, 2017), and even if they have not been highlighted in our screen, they might also regulate ERK dynamics. While providing the experimental throughput to perturb and analyze ERK dynamics at scale, optoFGFR, that lacks an ectodomain, evoked different ERK dynamics than endogenous RTKs such as FGFR and EGFR (Fig 3A and B compared to Fig 2F). These different ERK dynamics emerge likely because of receptor-level interactions that involve competition of bFGF for FGFR and heparan sulfate proteoglycan co-receptors (Kanodia et al, 2014;Blum et al, 2019) in the case of FGFR, or receptor endocytosis in the case of EGFR (Gerosa et al, 2020;Kiyatkin et al, 2020). Our combined modeling and experimental approach suggested that optoFGFR and EGFR share similar downstream MAPK network circuitries and NFBs (Fig 3C-G). OptoFGFR therefore provides a simplified system that allowed us to focus on intracellular feedback structures, without confounding receptor-level regulations. Our Bayesian inference modeling approach, that is parameter agnostic, could provide simple intuitions about the receptor-level and negative feedback structures that shape ERK dynamics in response to optoFGFR and EGFR inputs. However, even if we had access to many ERK dynamics phenotypes, our modeling approach did not allow us to explore more sophisticated MAPK network topologies such as the presence of two NFBs or multiple node isoforms. We interpreted our data using some of the feedback structures that have been previously experimentally documented and modeled but cannot formally exclude that the observed ERK dynamics emerge from different network structures. In the future, information about the different nodes and their dynamics might allow to further constrain the model topology and parameter space, and hopefully address this limitation.

Additional novel insights into regulation of ERK dynamics
Our optoFGFR and optoSOS screens provided new system-wide insights into the regulation of the MAPK network. Strikingly, the same perturbations induced different ERK dynamics phenotypes in the optoFGFR and optoSOS screens. This might occur because some regulators target the MAPK network at multiple levels, differently affecting ERK responses triggered with optoFGFR or optoSOS inputs. Additionally, as the two optogenetic systems are under the regulation of one versus two simultaneously occurring NFBs, they might have different sensitivities to perturbations, as discussed above.
With respect to the optoFGFR system, GRB2 KD led to a reduction of ERK amplitude (Fig 4D and E). GRB2 acts as the RTKproximal adaptor to activate SOS (Chardin et al, 1993;Belov & Mohammadi, 2012). As GRB2 operates at the start of the cascade, outside of most NFBs, heterogeneity in its expression levels might be less easily buffered out. PLCG1 KD increased damped oscillatory behavior (Fig 4E and F). Phospholipase Cɣ1 activates calcium signaling, which has itself been shown to regulate RAS/MAPK signaling in a calcium spike frequency-dependent manner (Cullen & Lockyer, 2002;Kupzig et al, 2005). Further investigation will be required to understand the significance of this crosstalk. RKIP KD resulted in higher ERK baseline and slower ERK adaptation poststimulation, without affecting ERK amplitude ( Fig 4D). RKIP (RAF kinase inhibitor protein) prevents MEK phosphorylation by CRAF (Yeung et al, 2000), suggesting that RKIP-dependent regulation is specifically involved in keeping basal ERK activity low. With respect to phosphatases, none of their perturbations led to a strong phenotype such as sustained ERK dynamics post-stimulation for example. The strongest phenotype was observed for PP2A KD that led to increased ERK amplitude, baseline, and slower adaptation (Figs 4D and EV3A). This might occur because the protein phosphatase 2A is an ubiquitous phosphatase that acts at multiple levels by dephosphorylating SHC1, MEK1, MEK2, ERK1 and ERK2, as well as a large number of other proteins (Junttila et al, 2008;Saraf et al, 2010). The observation that in optoFGFR-low PP2A KD cells, ERK dynamics displayed increased amplitude but still oscillated rather than exhibiting sustained behavior, suggests that NFBs might buffer against the loss of phosphatase regulation to some extent. Perturbation of the nuclear DUSPs, DUSP1,2,4, the atypical DUSP3 and most strongly the cytosolic DUSP6 (Patterson et al, 2009) led to higher ERK baseline, reduced adaptation, with only limited effects on amplitude (Figs 4D and EV3A). Consistently, DUSP6 has previously been proposed to pre-emptively dephosphorylate MAPKs to maintain low ERK activity baseline levels at resting state (Huang & Tan, 2012). Our results indicate that perturbation of single DUSPs might not be compensated by the others, suggesting that individual DUSPs might regulate specific substrates within the MAPK network. Except for DUSP6, KD of the different DUSPs did not significantly affect oscillatory ERK behavior in optoFGFR-low cells (Fig 4F), suggesting that they are not involved in the MAPK feedback circuitry that operates on timescales of minutes.
The optoSOS screen revealed stronger ERK amplitude phenotypes, especially for ERK2 and CRAF KD (Fig 5F versus Fig 4D). Unlike for optoFGFR input, RSK2 KD did not result in slower ERK adaptation, suggesting that ERK responses triggered by the optoSOS input are not regulated by the RSK2-mediated NFB. However, RSK2 KD led to a reduction of ERK amplitude, also observed to a lesser extent in response to optoFGFR input, suggesting a role of RSK2 in ERK amplitude regulation downstream of RAS. With respect to phosphatases, PP2A KD led to decreased amplitude, a different phenotype than in response to optoFGFR input. This might occur because of the broad specificity PP2A phosphatase, which might lead to different phospho-proteomes in response to optoSOS versus optoFGFR input. Similar phenomena might apply for most of the DUSPs.
The RSK2-mediated feedback can be targeted to potently inhibit oncogenic ErbB2 signaling Our data suggest that the RSK2-mediated NFB is important for MAPK signaling robustness downstream of our prototypic optoFGFR RTK (Fig 6). We found that the RSK2-mediated NFB likely also operates downstream of EGFR and oncogenic ErbB2 signaling in MCF10A cells (Fig 7). In response to EGF stimulation, or ErbB2 overexpression, a subset of RSK-inhibited cells displayed wider ERK pulses, suggesting that the RSK2 NFB is also involved in ERK adaptation in this system (Fig 7G cluster 3, Fig 7K cluster 1 and 2). Further, RSK inhibition led to a high heterogeneity of ERK dynamics within the cell population especially visible in the case of ErbB2 overexpressing cells (Fig 7J), which might result from the reduced ability of the MAPK network to cope with nodes expression noise in absence of the RSK2 NFB. In EGF-treated cells, combination of RSK and suboptimal MEK inhibition led to strong and homogeneous ERK inhibition (Fig 7E-G, cluster 5). In the ErbB2 overexpressing cells, combined RSK/MEK inhibition shifted most of the cell population to flat, low amplitude ERK dynamics, enabling to further inhibit a large number of cells when compared to suboptimal MEK inhibition only (Fig 7I-K, cluster 5). These results suggest that pharmacological inhibition of the RSK2-mediated NFB can be used to reduce MAPK signaling robustness, sensitizing the network to MEK perturbation. Such non-trivial drug combinations might allow for homogeneous inhibition of ERK dynamics in most of the cells in a population. This homogeneous inhibition might mitigate the emergence of drugtolerant persister cells from cell subpopulations that display residual ERK activity in response to inhibition of a single node. Our results imply that efficient pharmacological inhibition of the MAPK network requires precise understanding of its topology. The RSK2 NFB is an example of a druggable node that can be exploited to target MAPK signaling robustness.
Our scalable experimental pipeline provides new insight into the MAPK network wiring that produces ERK dynamics. However, our perturbation approach only highlighted very subtle ERK dynamics phenotypes, precluding a complete understanding of the MAPK network. We envision that this will require more precise knowledge about the dynamics of MAPK network nodes and their interactions in response to defined inputs and perturbations. Such data can now be produced at scale using optogenetic actuator/biosensor circuits as those we describe in this work. This information might allow for faithful parametrization of more complex models. With the increasing amount of optogenetic actuators and biosensors available, similar genetic circuits could also be designed to study the dynamics of other signaling pathways at scale.

Live imaging of ERK dynamics
NIH3T3 cells were seeded in 96-well 1.5 glass-bottom plates (Cellvis) coated with 10 μg/ml Fibronectin (Huber lab) using 1.5 × 10 3 cells/well and incubated for 24 h. MCF10A cells were seeded in 24-well 1.5 glass-bottom plates (Cellvis) coated with 5 μg/ ml Fibronectin (Huber lab) at 1 × 10 5 cells/well and incubated for 48 h. NIH3T3 cells were washed with PBS and incubated in starving medium for 4 h in the dark before starting the experiment. MCF10A cells were starved for 7 h before starting the experiments. In experiments involving drug perturbations, cells were incubated for 2 h (or 1 h in MCF10A experiments) with the inhibitor(s). Imaging was performed with an epifluorescence Eclipse Ti inverted fluorescence microscope (Nikon) using a Plan Apo air 20× (NA 0.8) objective. Nikon Perfect Focus System (PFS) was used to keep cells in focus throughout the experiment. Illumination was done with a SPECTRA X light engine (Lumencor) with the following filters (Chroma): mTurquoise2: 440 nm LED, 470lp, 69308 CFP/YFP/mCherry-ET, CFP 458-482; mCitrine: 508 nm LED, ET500/20×, 69308bs, ET535/ 30m; mRuby2 and mCherry: 555 nm LED, ET575/25×, 69008bs, 59022m, miRFP703: 640 nm LED, ET640/30×, 89100bs Sedat Quad, 84101m Quad. Images were acquired with an Andor Zyla 4.2 plus camera at a 16-bit depth. Image acquisition and optogenetic stimulation were controlled with the NIS-Element JOBS module. For NIH3T3 experiments, ERK-KTR-mRuby2 and H2B-miRFP703 were acquired at 1-min interval and 470-nm light inputs were delivered at specific frequencies and intensities (see below). MCF10A image acquisition was performed at 5-min time resolution. Growth factor stimulations were done by manually pipetting EGF and bFGF during the experiment. We used mCitrine intensity to quantify the expression level of the optogenetic constructs. However, as mCitrine excitation leads to optoFGFR or optoSOS activation, we acquired one frame with the ERK-KTR-mRuby2, the H2B-miRFP703, and the mCitrine-tagged optoFGFR or optoSOS only at the end of each NIH3T3 experiments. All experiments were carried on at 37°C with 5% CO 2 .

Optogenetic stimulation
Light stimulations were delivered with a 470 nm LED light source that was hardware-triggered by the camera to generate light pulses of reproducible duration. Light stimulations of defined intensity and duration were programmed to be automatically delivered at specific timepoints. To define the dose of light received by the cells, we measured the 470 nm light intensity at the focal plane using an optical power meter (X-Cite Power Meter, Lumen Dynamics Group) and converted this value to a power density as with FN = 18 mm. The obtained value was then multiplied by the duration of the pulse to obtain the dose of light received by the cells for each light pulse.
Light dose ðDÞ ¼ Light power density Â Pulse duration For stimulation of the optoFGFR cells, the 470-nm LED intensity was limited to a low dose by combining a ZET470/10× filter and a ND filter 5% (Chroma). Transient stimulations were done with a single pulse, while sustained stimulations were done with single pulses delivered every 2 min. For stimulation of the optoSOS cells, we used the 470 nm LED with a ET470/24× filter (no ND filter). Transient stimulations were done with 5 pulses repeated at 20-s intervals, while sustained stimulations were done using 5 pulses repeated at 20-s intervals, delivered every 2 min.

TIRF imaging of optoFGFR dynamics
Cells were seeded at a density of 1 × 10 3 per well in 96 well 1.5 glass bottom plates (Cellvis) coated with 10 μg/ml Fibronectin (Huber lab) and incubated for 24 h at 37°C with 5% CO 2 . Before imaging, cells were washed with PBS and incubated in starving medium for 4 h in the dark. Imaging was performed with an epifluorescence Eclipse Ti inverted fluorescence microscope (Nikon) using a CFI Apochromat TIRF 100× oil (NA 1.49). Images were acquired with an Andor Zyla 4.2 plus camera at a 16-bit depth. TIRF images were acquired with a 561 nm laser using a ET575/25 filter in front of the ZT488/561rpc (Chroma) to prevent nonspecific activation of the CRY2. MetaMorph software (Universal Imaging) was used for acquisition. TIRF images of the optoFGFR-mScarlet were acquired at a 20-s interval. Optogenetic stimulation was done using a 470 nm LED (SPECTRA X, Lumencor) (Appendix Fig S1B). All experiments were carried on at 37°C with 5% CO 2 .

Image processing pipeline
Nuclear segmentation was done in CellProfiler 3.0 (McQuin et al, 2018) using a threshold-based approach of the H2B channel. In the case of MCF10A cells, nuclear segmentation was preceded by prediction of nuclear probability using a random forest classifier based on different pixel features available in Ilastik software (Berg et al, 2019). To measure the ERK-KTR fluorescence in the cytosol, the nuclear mask was first expanded by 2 pixels to exclude the blurred edges of the nucleus. The new mask was then further expanded by 4 pixels in a threshold-based manner to obtain a "ring" area corresponding to the cytoplasmic ERK-KTR. ERK activity was obtained by calculating the ratio between the average cytosolic pixel intensity and the average nuclear pixel intensity. Single-cell tracking was done on nuclear centroids with MATLAB using μ-track 2.2.1 (Jaqaman et al, 2008). The final images containing the ERK-KTR-mRuby2, H2B-miRFP703, and the optoFGFR-mCitrine (or optoSOS-mCitrine) channels were processed using the same CellProfiler settings as the time-lapse images. Intensity of the mCitrine was extracted under the ERK-KTR cytoplasmic mask and used to classify cells into low or high expressors in a threshold-based manner. For optoFGFR-evoked ERK responses, the threshold was defined empirically to separate oscillatory and non-oscillatory ERK responses (low < −1.75 (log10 mCitrine intensity) < high). For optoSOS-evoked ERK responses, the threshold was defined empirically to separate cells with low or high ERK response amplitudes (low < −1.25 (log10 mCitrine intensity) < high). The same thresholds were kept across experiments to compare low and high expressors. The optoFGFR-mScarlet dimers/oligomers were segmented using the pixel classification module from Ilastik (Berg et al, 2019). OptoFGFR dimers, cell background and trafficking vesicles were manually annotated on images before and after the light stimulation. A probability map of the optoFGFR dimers classification was exported as TIFF for each frame. We then computed the mean of pixel intensities from the binarized mask obtained with Ilastik using Fiji (Appendix Fig S1C).

Quantification of ERK activity
We wrote a set of custom R scripts to automatically calculate the ERK-KTR cytoplasmic/nuclear ratio as a proxy for ERK activity for each single-cell, link single-cell ERK responses with the corresponding optoFGFR/optoSOS intensity value and export the corresponding ERK single-cell trajectories. For NIH3T3 data, outliers in ERK single-cell trajectories were removed using a clustering-based approach (https://github.com/pertzlab/Outlier_app). Trajectories with an ERK activity higher than 0.8 or lower than 0.2 before stimulation, above 1.6 during the whole experiment or displaying single time point spiking values were removed. For MCF10A data, trajectories with an ERK activity above 2 or shorter than 90% of the total experiment duration were removed. All the R codes used for further analysis are available as supplementary information (see Data availability section). Hierarchical clustering analysis of single-cell trajectories (Figs 2D,EV3F and G,and EV4B) was done using Time Course Inspector (Dobrzyń ski et al, 2019).

Modeling
The model for the EGF and light stimulated ERK cascade is a kinetic model, representing the EGF receptor, the inter-cellular proteins (RAS, RAF, MEK, ERK) as well as a negative feedback (NFB) from ERK to RAF and the inactivation of the EGF receptor in the form of endocytosis (Fig 3D). We explicitly modeled the ERK-KTR readout through nuclear and cytosolic KTR. The initial fraction of cytosolic KTR is estimated from the data through the parameter ktr init . The KTR readout Y(t) was taken to be the ratio of cytosolic KTR over nuclear KTR with additive Gaussian noise.
where the variance of the measurement noise σ 2 was estimated from the data. Appendix Table S1 shows all modeled species, their notation used for the equation, as well as the initial values. We assume that in the beginning of the experiment, all species are in the inactive form, reflecting the fact that the cells have been starved. The total concentrations of all species have been normalized to 1. The model equations are shown in Appendix Table S2. The phosphorylation events are modeled with Michaelis-Menten kinetics. The NFB is modeled through the modeling species NFB and its "active" version NFB*, which affects the dephosphorylation rate of RAF linearly. The activation, endocytosis, and recycling of the EGF receptor is modeled linearly. The model parameters are described in Appendix  Table S3. For the modeling of the two smaller models (without feedback ( Fig EV1E) or without endocytosis (Fig EV1H)), we set the corresponding parameters (k nfb and r 2,3 ) to zero.
For the parameter inference, we used a Nested Sampling algorithm as described in (Mikelson & Khammash, 2020). The inference was performed on the ETH High-performance Cluster Euler and was done using the parallel implementation on 48 cores. The algorithm was run for 24 h or until the algorithm stopped because the termination criterion Δ LFNS (see (Mikelson & Khammash, 2020) for details) was À∞. As prior distributions, we chose for all parameters noninformative log-uniform priors between 10 −5 and 10 5 , except for ktr init for which we chose a uniform prior on the interval [0, 1] and for σ for which we chose a log-uniform prior between 10 −5 and 1. Predictive distributions can be found on Figs 3F and G, and EV1F, G, I and J.

RNAi perturbation screen
We used Ingenuity Pathway Analysis (IPA, Qiagen) to select proteins directly interacting with ERK, MEK, RAF, RAS and FGFR, that are known to be expressed in NIH3T3 cells using a proteomics approach ( Table S4). We then imported this protein list in STRING (Jensen et al, 2009) to generate an interaction network with a minimum interaction score of 0.4. The final interactome was manually modified to display the protein names to facilitate the readout (Fig 4A). We targeted these selected proteins with RNA interference, using the siPOOL technology (one siPOOL containing a mix of 30 siRNAs targeting the same gene (Hannus et al, 2014), sequences available in the Data availability section). We arranged the siPOOLs in a 96 well plate format (in columns 2-5 and 8-11, one well per siPOOL) with the non-targeting siRNA (CTRL) and the positive control (mix of 5 nM siPOOL against ERK1 and 5 nM siPOOL against ERK2) placed alternately in columns 1, 6, 7, and 12. Cells were reverse transfected using RNAiMAX (Thermo Fisher, 13778150) following the recommended siPOOL transfection protocol (https:// sitoolsbiotech.com/protocols.php). OptoFGFR-expressing cells were transfected with 10 nM of siPOOL in a 96-well 1.5 glass-bottom plate (Cellvis) coated with 10 μg/ml Fibronectin (Huber Lab) at 0.3 × 10 3 cells/well density and incubated for 72 h at 37°C and 5% CO 2 . For the imaging, the 96-well plate was divided into 15 subexperiments, each sub-experiment consisting of a negative control well, a positive control well and 4 wells with different siPOOLs. We selected 2 FOVs per well and programmed the microscope to run the 15 experiments sequentially, acquiring the ERK-KTR-mRuby2 and the H2B-miRFP703 channels with a 1-min interval, stimulating the cells with sustained optoFGFR input (2-min intervals, D = 18 mJ/cm 2 ), and acquiring a final frame with ERK-KTR-mRuby2, H2B-miRFP703, and optoFGFR-mCitrine (Figs 4B and D-F, EV2C and D, EV3A-D, 6F-H and EV5B). For the optoSOS system, we limited the perturbation screen to targets acting below RAS (Figs 5F and G, and EV4D-G). Stimulations were done with sustained optoSOS input (5 repeated pulses at 2-min intervals, D = 0.6 J/ cm 2 ). For EGF experiments, cells were stimulated with 1 ng/ml EGF at t = 5 min (Fig EV3E-G).

Real-time qPCR
Cells were transfected with different concentrations of siPOOL in a 24-well plate at 5 × 10 3 cells/well density and incubated at 37°C with 5% CO 2 for 72 h before RNA isolation. Reverse transcription was done with the ProtoScript II reverse transcriptase kit (Bioconcept, M0368L). Real-time qPCRs were run using the MESA Green pPCR MasterMix Plus for SYBR Green assay (Eurogenetec, RT-SY2X-03+WOU) on the Rotor-Gen Q device (Qiagen). Each sample was tested in triplicate. Expression level of the gene of interest was calculated using the 2 −ΔΔCt method with GAPDH expression level as internal control (Fig EV2A). The following primers were used for the RT-qPCR (designed with the Real-time PCR (TaqMan) Primer and Probes Design Tool from GenScript).

Target
Forward sequences Reverse sequences Immunoblotting Cells were transfected with 10 nM siPOOL in 6 well plates at 6 × 10 4 cells/well density and incubated at 37°C with 5% CO 2 for 72 h. Cells were lysed in a buffer containing 10 mM

Time series feature extraction
We used custom scripts to extract features of ERK responses to transient optoFGFR input (Figs 2B and EV1A and B), sustained GF input ( Fig EV1C and D), and transient optoSOS input (Fig EV4A). The maximum peak (maxPeak) is the absolute value of the highest ERK activity in the trajectory. To estimate the full width at half maximum (FWHM), we first removed the baseline of the trajectories and increased their sampling frequency by a factor 30 with spline interpolation. On the resulting trajectory, we applied a "walk" procedure to quantify the FWHM. In this method, a pointer walks left and right (i.e., opposite and along the direction of time, respectively) from the maximum point of the trajectory. The pointer stops whenever the half maximum value is crossed. Both stops define a left and a right border, the time difference between these 2-border timepoints gives the FWHM. To avoid reporting aberrant FWHM values in cases where a peak cannot be clearly defined, we excluded FWHM calculation for trajectories where the fold change between the baseline (mean activity before stimulation) and the maximum value of the trajectory was below a threshold manually defined. ERKpostStim is the absolute value of ERK activity extracted 9 min after the last stimulation pulse to evaluate ERK adaptation. Statistical analysis ( Fig  EV1A and B) was done by comparing all conditions to the 20-min interval stimulation patterns with a Wilcoxon test using the FDR P-value correction (NS: non-significant, *< 0.01, **< 0.001, ***< 0.0001, and ****< 0.00001).
To evaluate ERK phenotypes under siRNA perturbations in response to sustained optoFGFR or optoSOS input (Figs 4D, EV3A, 5F and EV4F), we extracted the baseline (average ERK activity on 5 timepoints before stimulation), the maxPeak (maximum ERK activity within a 10-min time window following the start of the stimulation) and the ERKpostStim (ERK activity at a fixed timepoint post-stimulation (t optoFGFR = 42 min and t optoSOS = 40 min)) from 3 technical replicates. To avoid heterogeneity due to differences in optogenetic expression, we focused our analysis on cells with high optogenetic expression. The obtained baseline, maxPeak and ERKpostStim for each siRNA perturbation were z-scored to the nontargeting siRNA (CTRL). Non-significant results were manually set to gray. Statistical analysis was done by comparing each perturbation to the control with a Wilcoxon test using the FDR P-value correction (NS: non-significant, *< 0.05, **< 0.005, ***< 0.0005, ****< 0.00005).
For the comparison of both optogenetic systems (Figs 5E and EV4C), ERK baseline was obtained by averaging ERK activity on 5 timepoints before stimulation and ERK maxPeak was extracted within a 10-min time window following the start of the stimulation. Statistical analysis was done by comparing low and high expressing cells within and across optogenetic systems with a Wilcoxon test using the FDR P-value correction (NS: non-significant, *< 0.05, **< 0.005, ***< 0.0005, ****< 0.00005).
To quantify the efficiency of the three MAPK inhibitors on the reduction of ERK amplitudes under sustained high optoFGFR or optoSOS input (Fig 6), extraction of the maxPeak was limited by the fact that several concentrations led to a full suppression of ERK amplitudes. Therefore, we extracted ERK amplitudes at a fixed time point following the start of the stimulation (t fixed optoFGFR = 15 min, t fixed optoSOS = 10 min). The obtained ERK amplitudes were then plotted for each concentration for a fixed number of cells randomly selected (Figs 6B and F,and EV5D). To calculate the EC 50 of each drug, we normalized the data by setting the mean ERK responses of the non-treated condition to 1 and the mean ERK responses of the maximum concentration to 0. EC 50 then was calculated by fitting a Hill function to the mean ERK activity of each concentration (Figs 6C and G, and EV5E, Appendix Tables S5-S7). The heterogeneity of ERK amplitude at the fixed time point was evaluated by computing the normalized standard deviation of the extracted ERK activity per condition (Figs 6D and H,and EV5F).

Identification of ERK dynamics phenotypes using CODEX
To investigate ERK dynamics phenotypes to siRNA perturbations, we first trained a convolutional neural network (CNN) to classify input ERK trajectories into any of the siRNA-perturbed conditions (Fig EV3C). For this purpose, we used a CNN architecture composed of 4 1D-convolution layers with 20 kernels of size 5, followed by a convolution layer with 20 kernels of size 3 and one layer of 10 kernels of size 3. The responses are then pooled with global average pooling to generate a vector of 10 features that is passed to a (10,63) fully connected layer for classification. Each convolutional layer is followed by ReLU and batch normalization. The CNN was trained to minimize the cross-entropy loss, with L2 weight penalty of 1e −3 .
To identify siRNA treatments that induced a distinctive phenotype, we selected the 10 conditions for which the CNN classification precision was the highest on the validation set (Appendix Table S4, "CODEX accuracy"). To these 10 conditions, we also added the negative control (non-targeting siRNA (CTRL)). We trained a second CNN, with the same architecture and training parameters, but limited to recognizing the 11 selected treatments to obtain a clear embedding of these hits. With this new model, we extracted the features used for the classification of the trajectories (i.e., the input representation after the last convolution layer) and projected them with tSNE (Python's sklearn implementation, perplexity of 100, learning rate of 600 and 2,500 iterations) (Fig EV3D). We selected 10 prototype curves for each treatment by taking the trajectories for which the second CNN's classification confidence (i.e., the probability for the actual class of the inputs) were the highest in the validation set ( Fig 4E, "CODEX").
To visualize the ERK dynamics landscape in MCF10A WT cells and in MCF10A cells overexpressing ErbB2, we trained one CNN for each cell line. These CNNs were trained to recognize the drug treatment applied on cells, using single-cell ERK traces as input. The architecture of the CNNs is the same as described previously. The only difference lies in the number of outputs in the final fully connected layer, which were set to the number of drug treatments. Features used for the classification of the trajectories were then projected with tSNE (Figs 7D and H,and EV5J and K).
To identify clusters gathering similar ERK dynamics (Fig 7E, F, I and J), we clustered trajectories based on their CNN features using a partition around medoids (PAM). This iterative algorithm is similar to K-means clustering. PAM defines the cluster centers (i.e., the medoids) as the observed data points which minimize the median distances to all other points in its own cluster. This makes PAM more robust to outliers than K-means which uses the average coordinates of a cluster to define its center. Representative trajectories were obtained by taking the medoids of each cluster and their four closest neighbors (Fig 7G and K). Distances between points were defined with the Manhattan distance between the scaled CNN features (zero mean and unit variance). We manually verified that these clusters captured an actual trend by visualizing trajectories in each cluster with the interactive CODEX application.

Peak detection and classification of oscillatory trajectories
The number of ERK activity peaks was calculated with a custom algorithm that detects local maxima in time series. First, we applied a short median filter to smoothen the data with a window width of 3 timepoints. Then, we ran a long median filter to estimate the longterm bias with a window width of 15 timepoints. This bias was then subtracted from the smoothed time series, and we only kept the positive values. If no point in this processed trajectory was exceeding a manual threshold of 0.075, all variations were considered as noise and no peak was extracted from the trajectory. The remaining trajectories were then rescaled to [0,1]. Finally, peaks were detected as points that exceeded a threshold which was manually set to 0.1. Peaks that were found before the first stimulation or after the last stimulation were filtered out.
The classification of trajectories into oscillatory and nonoscillatory behaviors was performed after the peak detection step. Cells were called oscillatory if at least 3 peaks were detected with the peak detection procedure (Fig 4F). Statistical analysis was done using a pairwise t-test comparing each perturbation to the control for high and low levels of optoFGFR independently, with FDR P-value correction (*< 0.05, **< 0.005, ***< 0.0005, ****< 0.00005).